tests/testthat/profiling for y_bootstrap_komplettfunktion.R

In nikosbosse/SAE: This Package implements a Small Area Estimation approach from Elbers. et al blabla

load("C:/Users/nikos/Desktop/SAE_proj_data.RData")

library(microbenchmark)
library(parallel)
library(mvtnorm)


# Alternative 1 ###############################################################################
# 2 schritte
bootstrap.y <- function(model1, model_fit1, censusdata1, n_boot1, n_obs){
  vars <- all.vars(model1)[-1] # wir können ja nochmal gucken, welche Variante richtig ist
  x_census <- as.matrix(cbind(rep(1, times = n_obs), subset(censusdata1, select = vars)))
  tx_census <- t(x_census)

  cov_coefficients <- vcov(summary(model_fit1))
  var_y <- rep(NA, n_obs)
  for(i in 1:n_obs){
    var_y[i] <- tx_census[,i] %*% cov_coefficients %*% x_census[i,]
  }
  sd_y <- sqrt(var_y)

  xbeta <- predict(model_fit1, newdata = censusdata1)

  y_bootstrap <- matrix(NA, nrow = n_obs, ncol = n_boot1)
  for (i in 1:n_obs){

    y_bootstrap[i,] <- rnorm(n = n_boot1, mean = xbeta[i], sd = sd_y[i])

  }

  return(y_bootstrap)

}

# Alternative 2 ###############################################################################
# ein schritt
bootstrap.y2 <- function(model1, model_fit1, censusdata1, n_boot1, n_obs){

  beta_bootstrap <- t(MASS::mvrnorm(n = n_boot1,
                                    mu = model_fit1$coefficients,
                                    Sigma = vcov(summary(model_fit1))))

  x_census <- as.matrix(cbind(rep(1, times = n_obs),
                              subset(censusdata1, select = all.vars(model1)[-1])))

  return(x_census %*% beta_bootstrap)
}

# Alternative 2.2 ###############################################################################
# ein schritt, Multiplikation mit C++
bootstrap.y2.2 <- function(model1, model_fit1, censusdata1, n_boot1, n_obs){

  beta_bootstrap <- t(MASS::mvrnorm(n = n_boot1,
                                    mu = model_fit1$coefficients,
                                    Sigma = vcov(summary(model_fit1))))

  x_census <- as.matrix(cbind(rep(1, times = n_obs),
                              subset(censusdata1, select = all.vars(model1)[-1])))

  return(matrixproductC(x_census, beta_bootstrap))
}



mbm <- microbenchmark::microbenchmark(
  alt1 = bootstrap.y(model1 = model, model_fit1 = inference_survey$model_fit_surv, censusdata1 = censusdata,
                     n_boot1 = 250, n_obs = n_obs_census),
  alt2 = bootstrap.y2(model1 = model, model_fit1 = inference_survey$model_fit_surv, censusdata1 = censusdata,
                     n_boot1 = 250, n_obs = n_obs_census),
  times = 10);mbm


# Alternative 2 is a lot faster

mbm <- microbenchmark::microbenchmark(
  alt2 = bootstrap.y2(model1 = model, model_fit1 = inference_survey$model_fit_surv, censusdata1 = censusdata,
                     n_boot1 = 250, n_obs = n_obs_census),
  alt2.2 = bootstrap.y2.2(model1 = model, model_fit1 = inference_survey$model_fit_surv, censusdata1 = censusdata,
                      n_boot1 = 250, n_obs = n_obs_census),
  times = 2);mbm

# seems like Alt 2 is still faster